#ifndef __OPERATOR_H__
#define __OPERATOR_H__

#include "Equation.h"

template <int DIM>
class RestrictionOperator
{
protected:
	int n;

public:
	RestrictionOperator(int _n);
	virtual ~RestrictionOperator() = default;
	virtual const Vector operator*(const Vector& x) = 0;
};

template <int DIM>
RestrictionOperator<DIM>::RestrictionOperator(int _n) : n(_n) {}



template <int DIM>
class FullWeighting : public RestrictionOperator<DIM>
{
public:
	FullWeighting(int _n);
	~FullWeighting() = default;
	virtual const Vector operator*(const Vector& x);
};

template <int DIM>
FullWeighting<DIM>::FullWeighting(int _n) : RestrictionOperator<DIM>(_n) {}

template <int DIM>
const Vector FullWeighting<DIM>::operator*(const Vector& x)
{
	Vector y;
	int n = RestrictionOperator<DIM>::n;
	if (DIM == 1) {
		if (x.size() != n - 1)
			throw std::length_error("FullWeighting: sizeof(x) != n-1");
		y = Vector(n / 2 - 1);
		for (int i = 0; i <= n / 2 - 2; ++i)
			y[i] = (x[2 * i] + 2 * x[2 * i + 1] + x[2 * i + 2]) / 4.0;
	}
	else if (DIM == 2) {
		if (x.size() != (n - 1) * (n - 1))
			throw std::length_error("FullWeighting: sizeof(x) != (n-1)^2");
		y = Vector((n / 2 - 1) * (n / 2 - 1));
		int k0, k1, k2, k3, k4, k5, k6, k7, k8;
		for (int i = 2; i <= n - 2; i += 2) {
			for (int j = 2; j <= n - 2; j += 2) {
				k0 = (n - 1) * (j - 2) + (i - 2);
				k1 = (n - 1) * (j - 2) + (i - 1);
				k2 = (n - 1) * (j - 2) + i;
				k3 = (n - 1) * (j - 1) + (i - 2);
				k4 = (n - 1) * (j - 1) + (i - 1);
				k5 = (n - 1) * (j - 1) + i;
				k6 = (n - 1) * j + (i - 2);
				k7 = (n - 1) * j + (i - 1);
				k8 = (n - 1) * j + i;
				y[(n / 2 - 1) * (j / 2 - 1) + i / 2 - 1] =
					(x[k0] + x[k1] + x[k2] + x[k3] + 8.0 * x[k4] + x[k5] + x[k6] + x[k7] + x[k8]) / 16.0;
			}
		}
	}

	return y;
}



template <int DIM>
class Injection : public RestrictionOperator<DIM>
{
public:
	Injection(int _n);
	~Injection() = default;
	virtual const Vector operator*(const Vector& x);
};

template <int DIM>
Injection<DIM>::Injection(int _n) : RestrictionOperator<DIM>(_n) {}

template <int DIM>
const Vector Injection<DIM>::operator*(const Vector& x)
{
	Vector y;
	int n = RestrictionOperator<DIM>::n;
	if (DIM == 1) {
		if (x.size() != n - 1)
			throw std::length_error("Injection: sizeof(x) != n-1");
		y = Vector(n / 2 - 1);
		for (int i = 0; i <= n / 2 - 2; ++i)
			y[i] = x[2 * i + 1];
	}
	else if (DIM == 2) {
		if (x.size() != (n - 1) * (n - 1))
			throw std::length_error("FullWeighting: sizeof(x) != (n-1)^2");
		y = Vector((n / 2 - 1) * (n / 2 - 1));
		for (int i = 2; i <= n - 2; i += 2) {
			for (int j = 2; j <= n - 2; j += 2) {
				y[(n / 2 - 1) * (j / 2 - 1) + i / 2 - 1] = x[(n - 1) * (j - 1) + (i - 1)];
			}
		}
	}

	return y;
}



template <int DIM>
class InterpolationOperator
{
protected:
	int n;

public:
	InterpolationOperator(int _n);
	virtual ~InterpolationOperator() = default;
	virtual const Vector operator*(const Vector& x) = 0;
};

template <int DIM>
InterpolationOperator<DIM>::InterpolationOperator(int _n) : n(_n) {}



template <int DIM>
class Linear : public InterpolationOperator<DIM>
{
public:
	Linear(int _n);
	~Linear() = default;
	virtual const Vector operator*(const Vector& x);
};

template <int DIM>
Linear<DIM>::Linear(int _n) : InterpolationOperator<DIM>(_n) {}

template <int DIM>
const Vector Linear<DIM>::operator*(const Vector& x)
{
	Vector y;
	int n = InterpolationOperator<DIM>::n;
	if (DIM == 1) {
		if (x.size() != n / 2 - 1)
			throw std::length_error("sizeof(x) != n/2 - 1");
		y = Vector(n - 1);
		y[0] = x[0] / 2.0; y[n - 2] = x[n / 2 - 2] / 2.0;
		y[1] = x[0];
		for (int i = 1; i <= n / 2 - 2; ++i) {
			y[2 * i + 1] = x[i];
			y[2 * i] = (x[i] + x[i - 1]) / 2.0;
		}
	}
	else if (DIM == 2) {
		if (x.size() != (n / 2 - 1) * (n / 2 - 1))
			throw std::length_error("sizeof(x) != (n/2 - 1)^2");
		y = Vector((n - 1) * (n - 1));
		double value;
		for (int i = 2; i <= n - 2; i += 2) {
			for (int j = 2; j <= n - 2; j += 2) {
				value = x[(n / 2 - 1) * (j / 2 - 1) + (i / 2 - 1)];
				y[(n - 1) * (j - 1) + (i - 1)] += value;
				y[(n - 1) * (j - 2) + (i - 1)] += 0.5 * value;
				y[(n - 1) * (j - 1) + i] += 0.5 * value;
				y[(n - 1) * j + (i - 1)] += 0.5 * value;
				y[(n - 1) * (j - 1) + (i - 2)] += 0.5 * value;
				y[(n - 1) * (j - 2) + (i - 2)] += 0.25 * value;
				y[(n - 1) * (j - 2) + i] += 0.25 * value;
				y[(n - 1) * j + i] += 0.25 * value;
				y[(n - 1) * j + (i - 2)] += 0.25 * value;
			}
		}
	}

	return y;
}



template <int DIM>
class Quadratic : public InterpolationOperator<DIM>
{
public:
	Quadratic(int _n);
	~Quadratic() = default;
	virtual const Vector operator*(const Vector& x);
};

template <int DIM>
Quadratic<DIM>::Quadratic(int _n) : InterpolationOperator<DIM>(_n) {}

template <int DIM>
const Vector Quadratic<DIM>::operator*(const Vector& x)
{
	Vector y;
	int n = InterpolationOperator<DIM>::n;
	if (DIM == 1) {
		if (x.size() != n / 2 - 1)
			throw std::length_error("sizeof(x) != n/2 - 1");
		y = Vector(n - 1);
		if (n == 4) {
			y[0] = y[2] = 3.0 * x[0] / 4.0;
			y[1] = x[0];
		}
		else {
			y[0] = 3.0 * x[0] / 4.0 - x[1] / 8.0;
			y[2] = 3.0 * x[0] / 4.0 + 3.0 * x[1] / 8.0;
			y[n - 2] = -x[n / 2 - 3] / 8.0 + 3.0 * x[n / 2 - 2] / 4.0;
			y[n - 3] = x[n / 2 - 2];
			y[n - 5] = x[n / 2 - 3];
			for (int i = 0; i <= n / 2 - 4; ++i) {
				y[2 * i + 1] = x[i];
				y[2 * i + 4] = -x[i] / 8.0 + 3.0 * x[i + 1] / 4.0 + 3.0 * x[i + 2] / 8.0;
			}
		}
	}
	else if (DIM == 2) {
		if (x.size() != (n / 2 - 1) * (n / 2 - 1))
			throw std::length_error("sizeof(x) != (n/2 - 1)^2");
		y = Vector((n - 1) * (n - 1));
		double value;
		for (int i = 2; i <= n - 2; i += 2) {
			for (int j = 2; j <= n - 2; j += 2) {
				value = x[(n / 2 - 1) * (j / 2 - 1) + (i / 2 - 1)];
				y[(n - 1) * (j - 1) + (i - 1)] += value;
				y[(n - 1) * (j - 2) + (i - 2)] += 0.25 * value;
				y[(n - 1) * (j - 2) + i] += 0.25 * value;
				y[(n - 1) * j + i] += 0.25 * value;
				y[(n - 1) * j + (i - 2)] += 0.25 * value;
				y[(n - 1) * (j - 2) + (i - 1)] += 0.5625 * value;
				y[(n - 1) * (j - 1) + i] += 0.5625 * value;
				y[(n - 1) * j + (i - 1)] += 0.5625 * value;
				y[(n - 1) * (j - 1) + (i - 2)] += 0.5625 * value;
				if (j - 3 >= 1)
					y[(n - 1) * (j - 4) + (i - 1)] -= 0.0625 * value;
				if (i + 3 <= n - 1)
					y[(n - 1) * (j - 1) + (i + 2)] -= 0.0625 * value;
				if (j + 3 <= n - 1)
					y[(n - 1) * (j + 2) + (i - 1)] -= 0.0625 * value;
				if (i - 3 >= 1)
					y[(n - 1) * (j - 1) + (i - 4)] -= 0.0625 * value;
			}
		}
	}

	return y;
}
#endif